arXiv:l502.05519vl [cond-mat.soft] 19 Feb 2015 


CREATED USING THE RSC LaTeX PCCP ARTICLE TEMPLATE - SEE www.rsc.org/electronicfiles FOR DETAILS 

ARTICLE TYPE www.rsc.org/xxxxxx | XXXXXXXX 


Collective waves in dense and confined microfluidic droplet arrays^ 

Ulf D. Schiller,* *^ Jean-Baptiste Fleury,^ Ralf Seemann,^^ and Gerhard Gompper"^ 

Received Xth XXXXXXXXXX 20XX, Accepted Xth XXXXXXXXX 20XX 
First published on the web Xth XXXXXXXXXX 200X 
DOI: 10.1039/b000000x 


Excitation mechanisms for collective waves in confined dense one-dimensional microfiuidic droplet arrays are investigated by 
experiments and computer simulations. We demonstrate that distinct modes can be excited by creating specific ‘defect’ patterns 
in flowing droplet trains. Excited longitudinal modes exhibit a short-lived cascade of pairs of laterally displacing droplets. 
Transversely excited modes obey the dispersion relation of microfiuidic phonons and induce a coupling between longitudinal 
and transverse modes, whose origin is the hydrodynamic interaction of the droplets with the confining walls. Moreover, we 
investigate the long-time behaviour of the oscillations and discuss possible mechanisms for the onset of instabilities. Our findings 
demonstrate that the collective dynamics of microfiuidic droplet ensembles can be studied particularly well in dense and confined 
systems. Experimentally, the ability to control microfiuidic droplets may allow to modulate the refractive index of optofiuidic 
crystals which is a promising approach for the production of dynamically programmable metamaterials. 


1 Introduction 

Microfiuidic devices have become an important tool in chem¬ 
istry and biology, where they are increasingly used, for ex¬ 
ample, in analytic essays,® micro-reactions® or flow cytom¬ 
etry.® These applications typically involve manipulation and 
control of immersed objects, such as droplets, vesicles or 
cells,® that interact hydrodynamically through the flow per¬ 
turbations of the surrounding fluid. A detailed understand¬ 
ing of the correlated motion induced by long-ranged hydro- 
dynamic interactions in microfiuidic devices is therefore es¬ 
sential for efficient control of the flow of micro-particles. In 
recent years, a number of studies brought out that microflu¬ 
idic droplet systems are especially well suited to steer their 
dynamics by modifying particle properties and/or device ge¬ 
ometry.®® Consequently, microfiuidic droplets have become 
both a test-bed and a model system to study collective be¬ 
haviour and self-organisation in non-equilibrium many-body 
systems.® Typically, a pressure-driven flow is imposed such 
that the system is out of equilibrium, and at low Reynolds 
number viscous dissipation dominates over inertia. A theo¬ 
retical description of such driven dissipative systems remains 
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challenging, and thus experiments and computer simulations 
are basic tools to study the dynamics of microfiuidic droplet 
ensembles. 

When microfiuidic droplets are confined between two par¬ 
allel plates, the geometry is effectively two-dimensional (2D) 
and the scattered flow has a characteristic dipolar form.® In 
this case, the hydrodynamic interactions are marginally long- 
ranged, i.e., the decay exponent is equal to the dimensionality 
of the system.® In contrast to quasi-ID geometries, where the 
hydrodynamic interactions are strongly screened, the dipolar 
interactions in quasi-2D geometries lead to complex collec¬ 
tive phenomena.® Dipolar flow fields are also characteristic 
for some types of self-propelled particles, such as droplets 
driven by Marangoni flows or by chemical reactions on their 
surface.® Some progress has been made in understanding the 
dynamics of rigid and deformable particles and their hydro- 
dynamic coupling in 2D pressure-driven flow. Pairs of rigid 
particles in Poiseuille flow were shown to follow either bound 
or unbound trajectories, depending on the relative position 
of the particles, their absolute position in a channel, and the 
strength of confinement.® Linear arrays of rigid spheres and 
deformable drops aligned in the flow direction undergo a pair¬ 
ing instability.® While arrays of spherical particles are also 
unstable to lateral perturbations, droplet arrays are stabilised 
by quadrupolar interactions due to deformation.® Asymmet¬ 
ric particles align with the flow due to self-interactions, and 
migrate to the centreline of the confining channel.®® For 
highly asymmetric particles, the time-scales for alignment and 
focusing separate due to the distinct hydrodynamic mecha¬ 
nisms involved. The focusing of asymmetric particles re¬ 
sembles a damped harmonic oscillator, whereas symmetric 
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particles oscillate between side-walls. ESI in 2D ensembles of 
droplets, the dipolar hydrodynamic interactions give rise to 
sound and shock waves that are superposed on droplet diffu¬ 
sion. The waves are due to a density-vel ocity coupling and 
can be described by a ID Burgers equation .E^E2] 

A particular feature arises in regular arrays of droplets, so 
called microfluidic crystals, where the flowing droplets have 
a spatial order with a well-defined spacing. These crystals 
can exhibit collective oscill ation s with a dispersion relation 
akin to solid state phonons.These microfiuidic phonon 
modes are neither growing nor decaying, and are thus a re¬ 
alisation of marginally stable oscillatory modes in a dissipa¬ 
tive system made possible by the imposed symmetry-breaking 
fiow. Practically, however, the possibili ty to observe these 
modes is limited by non-linear instabilitiesESEU strong 

dependence on initial conditions. Only recently, an experi¬ 
mental technique was proposed to systematically excite mi¬ 
crofiuidic phonons, and the observed modes revealed a cou¬ 
pling mechanism, induced by lateral confinement, between 
longitudinal and transverse modes that was confirmed by com¬ 
puter simulations .23 The ability to control the dynamic prop¬ 
erties by tuning the fiow characteristics opens interesting per¬ 
spectives regarding dynamically programmable metamaterials 
which could be produced by modulating the refractive index 
of droplet crystals.ESES 

Here we investigate collective modes in dense microfiuidic 
crystals under confinement both experimentally and by com¬ 
puter simulations. We show that distinct oscillatory behaviour 
can be systematically excited by varying the initial conditions 
through the introduction of specific ‘defect’ patterns. The ob¬ 
served modes are analysed and characterised, and reveal sev¬ 
eral interesting dynamic features, such as cascades of later¬ 
ally offset pairs and mode coupling. The results from exper¬ 
iments and computer simulations agree quantitatively. The 
long-time behaviour is investigated in computer simulations 
and used to identify possible instabilities and their underlying 
mechanisms. Our approach demonstrates the rich dynamics 
that emerges from hydrodynamic interactions in confined mi¬ 
crofiuidic droplet ensembles. The results show good agree¬ 
ment with a linearised far-field theoryESl even in the dense 
droplet regime. This makes it very promising to apply the 
techniques to other crowded microfiuidic systems, such as 
self-propelled particles. E3 

The remainder of the article is organised as follows: In sec- 
tion[^ we review the hydrodynamics of quasi-2D systems and 
the linearised far-field theory for microfiuidic phonons. Sec¬ 
tion describes the experimental techniques and the simula¬ 
tion approach we used to study microfiuidic droplet systems. 
In section]^ we present excitation mechanisms for collective 
waves and analyse the observed oscillations and instabilities. 
A concluding discussion is given in section 



Fig. 1 Schematic illustration of the quasi-2D flow geometry. 
Flattened droplets experience a friction at the top and bottom plates 
which counteracts the hydrodynamic drag. When the droplets move 
relative to the imposed flow, they act as a mass dipole with a sink at 
their leading edge and a source at their trailing edge. 


2 Microdroplet trains in quasi-2D flow 


We consider droplets that are confined between two parallel 
plates and thus move in a quasi-2D geometry, cf. figure 
The fiow and the hydrodynamic interactions in this geometry 
differ qualitatively from the bulk case due to momentum ab¬ 
sorption at the confining plates which leads to screening of the 
far-field. The fluid fiow satisfies no-slip boundary conditions 
on the channel walls, and since the height H of the channel is 
small compared to the lateral width W, the velocity gradient 
in the z-direction is much larger than in the planar directions. 
In the Darcy approximation ^ ^ dy, the solution 

of the Stokes equation has a quasi-2D Hele-Shaw formEIl 
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( 1 ) 


At low Reynolds number, the fiow is incompressible and 
Eq.Q can be written as a Laplace equation 

v20(x,j)=O (2) 

u(x,j) = V(/»(x,j) (3) 


where the effective potential is defined through the pressure 
(j) = —H^p/IIt] with 7] the dynamic viscosity of the fluid. 

We b riefly review here the theoretical description presented 
in Ref s.ESEI When a droplet is moving through the fluid with a 
velocity 5u = u°° — uj relative to the externally imposed fiow 
u°° = u°°x, it acts as a momentum monopole whose flux scales 
as 5u^. However, due to the absorption of momentum at the 
top and bottom plates the flux is not conserved. The absorbed 
flux scales as 5u/h, therefore the fiow field of the momen¬ 
tum monopole dr5u oc —5u/h decays exponentially .E3 Thus, 
unlike in the bulk case, the leading contribution is the mass 
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dipole created by the droplet. The moving droplet pushes fluid 
out at the upstream edge and draws fluid in at the downstream 
edge, giving rise to a characteristic dipolar flow fleld. For¬ 
mally, the flow perturbation at a distance r from a droplet is 
obtained by solving the Laplace equation with boundary con¬ 
ditions of zero mass flux (zero perpendicular velocity) through 
the droplet interface. EHl This gives the dipolar potential and 
scattered velocity fleld around a droplet of radius R 


and then rescaling by a compre ssibil ity factor C to account for 
the flnite size R of the droplets .CSIlll The result isEl 


7tR^5u 
+ coth 


,, ^ kr^ou f , r ^ 
w</W = C-^^|coth[—(z-%)J 


( 10 ) 


where 
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The 2D potential flow can also be described by a complex po¬ 
tential w(v + /y) = 0(v,y) + i\i/{x,y) where the imaginary part 
is the stream function i/A(x,y). The flow velocity is then given 
by Ux — iuy = dw/dz with z = x -1- iy. For an imposed flow in 
the v-direction the complex dipolar potential is then 


The equation of motion for a conflned droplet can be obtained 
as above and keeps the form = ldu°° if the mobility is re¬ 
placed by 

M=C-(c-^ + 0 . (12) 


Wd{z)=R^du^. ( 6 ) 

A droplet moving in the imposed flow experiences a hydro- 
dynamic drag that can be written as 

^h= + (7) 

where the drag coefficient ^ = 24717]R^/H is introduced. The 
second term arises from the self-interaction of the droplet with 
its dipole.nSI 

If the size R of the droplets exceeds the channel height H, 
they are flattened and experience a friction with the top and 
bottom plates which can be modelled as 


Ff=-Cuj. (8) 

Since inertial effects can be neglected at low Reynolds num¬ 
ber, we can use force balance Fh+Ff = 0 to obtain the equation 
of motion for the droplet 


Urf = ^u“ = Q + 0 «“> (9) 

where fl = Ud I u° is the mobility of the droplet in the imposed 
flow. 

In the presence of lateral side-walls, i.e., in a microfluidic 
channel, additional boundary conditions have to be satisfled. 
The simple dipole potential ([^ has a non-vanishing flux at the 
side-walls whi ch can be eliminated by placing image dipoles 
inside the wall.l^^lSl These dipoles form an inflnite array per¬ 
pendicular to the flow direction. The flow potential of a single 
droplet is obtained by summing over the inflnite dipole array. 


In an ensemble of droplets, the solution of the Laplace equa¬ 
tion is considerably more complicated because the boundary 
conditions have to be satisfled additionally on all droplet sur¬ 
faces. Although this is in principle possible using the method 
of images, the large number of reflections that arise makes it 
unreasonably intricate in practice. One therefore resorts to the 
leading-order approximation where the drag force is given by 
a superposition of the flow flelds created by the other droplets. 
For the ^-th droplet in an ensemble, the equation of motion 
thus is^SI 




jl \ U 
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dwd 

dz 


(13) 


This approximation is valid if the inter-droplet distance is 
larger than the droplet size rj — Vn ^ R, and we will see be¬ 
low that the predictions based on Eq. work well even for 
dense droplet trains. 

For a regular train of droplets flowing with an ‘equilib¬ 
rium’ spacing a in the centre of the channel, the displacements 
Szn =Zn—7ia are assumed to be small and the derivative of the 


potential ( p^ can be exp anded . To first order, the equations of 
motion are then given byllSE3 
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with prefactor B = — Ud){7r^R/W^)t3.n{7rR/W). These 

equations describe waves, and by plugging in a plane-wave 
solution, we arrive at the dispersion relations 


collision steps. During the streaming step the particle moves 
ballistically, 

r,-= r, + /zv,-, (16) 
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3 Experimental and computational droplet mi¬ 
crofluidics 

Having discussed the flow fields and the forces acting on a 
single droplet, we discuss the experimental realisation as well 
as computer simulations of flowing droplets in a straight mi- 
crofiuidic channel. 


3.1 Microchip fabrication and droplet production 

Microfiuidic devices were fabricated in Sylgard 1 84 (D ow 
Corning) using standard soft lithographic protocols,EEIl and 
the flow rates were volume-controlled by syringe pumps. 
Mono-disperse water droplets were generated in n-hexadecane 
(p = 773kg/m^, T] = 3mPas) with 2 wt% of the surfactant 
Span 80 using a step geometry,EEIl. The microchannel has 
uniform height and width of ^ x IT 120 pm x 210 pm. Typ¬ 
ical flow velocities are ~ 250 pm/s for the droplet, and 
Uoii ~ 500 pm/s for the continuous oil phase. The correspond¬ 
ing Reynolds and Peclet number are Re = pu^nR/T] 10“^ 
and Pe = UquR/D ^10^, respectively. 

3.2 Simulation approach: Multi-particle collision dy¬ 
namics 

In order to investigate the origin of our experimental obser¬ 
vations and the validity of the approximations in the theoreti¬ 
cal description, we conduct computer simulations using multi¬ 
particle collision dynamics (MPC). MPC is a mesoscopic sim¬ 
ulation method that is capable of reproducing the full hydro¬ 
dynamics of a fiuid.ISlESI 3 assump¬ 

tion of a specific flow perturbation, it is well suited to test the 
accuracy of the semi-analytical theory based on dipolar fiow 
fields to describe the droplet interactions in a dense and con¬ 
fined system. The fiuid is modelled explicitly by idealised 
point-like particles of mass m. The fiuid dynamics emerges 
from local mass, momentum and energy conservation in the 
particle ensemble, whose equation of state is that of an ideal 
gas. The update of particle positions and momenta mimics the 
underlying kinetics and is split into successive streaming and 


where h is the time interval between collisions. In the colli¬ 
sion step, the particles are sorted into cubic collision cells of 
size Ax. In each cell, the particles then exchange momentum 
while the momentum of the collision cell is conserved. Vari¬ 
ous collision rules have been proposed in the literature and in 
this work, we employ a collision rule that also conserves an¬ 
gular momentum of the cell.El The collisions are augmented 
with an Anderson-like thermostat to control the temperature. 
The overall update of particle velocities is 


y;=yc + yr-I^ 




+mn 


-1 


I 


rj,c X (y; 


“)] Vc, 


(17) 


where vc is the centre of mass velocity of the collision cell 
containing Nc particles, □ is the moment of inertia tensor of 
the particles, c = r/ — rc is the relative particle position, and 
v-^^ is a random velocity drawn from a Maxwell-Boltzmann 
distribution. This collision operator is denoted as MPC-AT-F (2 
in the nomenclature of Ref.^. In addition, the cell grid is 
shifted randomly before each collision step to restore Galilean 
invariance of the system.l^ The dynamic viscosity r] of the 
MPC-AT-F (2 fiuid for large number density n (particles per cell) 
is then given by 


nksTh / n 1\ m{n — ll5) 

^ Ax^ —(J + 2)/4 2y^ 24Ax^~^h 


(18) 


where ksT is the imposed temperature and J = 2 is the dimen¬ 
sionality of the system. 

Since the droplets hardly deform in the experiment, we 
model them as rigid discs of radius R that are coupled to 
the fiuid by a no-slip boundary condition, i.e., v' = —v + 2vb 
where Vb is the boundary velocity. It is to be noted that this is 
effectively a different boundary condition than the one used in 
deriving Eq. Q, however, we have found in practice that this 
can be accounted for by the calibration procedure described 
below and does not lead to a relevant difference in the mea¬ 
surements. To apply the collision rule in the cells that are 
partly or fully occupied by the rigid discs, the correspond¬ 
ing volume is filled with virtual particles that are distributed 
randomly within a layer of width ^/2a and whose velocities 
are distributed according to a Maxwell-Boltzmann distribution 
around the boundary velocity Vb-^ The momentum change 
of the fiuid particles during streaming and collisions is accu¬ 
mulated and leads to the boundary force Fb that moves the 

discs.®! 
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On very short distances, hydrodynamics are not resolved 
and hence we add steric droplet-droplet and droplet-wall in¬ 
teractions by means of Weeks-Chandler-Anderson (WCA) po¬ 
tentials 


VwcA{r) = 4e 



+ e 


0 < r < 26ro, 


(19) 

where r = rij — 2R for droplet-droplet interactions, and r = 
Vij — R for droplet-wall interactions. Vij is the distance of the 
droplet centres, or the distance of the droplet centre from the 
wall surface, respectively. To account for the friction at the top 
and bottom of the Hele-Shaw cell, cf. Eq. ([^, we apply a fric¬ 
tion force Ffriction = — 7 Ud to the discs where Ud is the velocity 
of the droplet relative to the microchannel. The overall flow 
is driven by an external force g corresponding to a constant 
pressure gradient across the channel. 

The parameters of the MFC simulations are as follows. 
The size of a collision cell is Ax = /?/5 and the time step is 
h = 0.005 T, where the time scale is T = {ksT /m)“^/^Ax, and 
m is the mass of one MFC particle. The fluid density, the driv¬ 
ing force and the friction are p = 40m/Av^, g = O.lmAx/T, 
and 7 = 2 • lO^m/T. These parameters correspond to a fluid 
viscosity 7] = 321.77m/T. The mass of the droplets is given 
by M = pTlR^ ^ 3140m. The parameters for the WCA poten¬ 
tial are £ = ksT and tq = Ax. We varied the channel width 
W, the droplet spacing a, and the initial conflguration includ¬ 
ing the initial wavelength A for simulating droplet trains in a 
channel. 

In order to compare the simulation results to the exper¬ 
iments quantitatively, we determined the value K = Udju^w 
from independent simulation runs with a single droplet un¬ 
der the same confinement. For a channel width of W = 3R 
we obtained a value ofK^ 0.62 which is on the order of 
the experimental parameters. The oil velocity is in the range 
Uoi\ « 10“^Ax:/t, such that the typical Reynolds and Feclet 
number of the simulations are Re ~ 10“^ and Pe ~ 10^, re¬ 
spectively. Note that the Feclet number is significantly lower 
than in the experiments. On the one hand, this leads to more 
pronounced fluctuations, but on the other hand, it allows us to 
observe in the simulations the onset of instabilities on acces¬ 
sible time scales, cf. section [4^ 


4 Controlled excitation and analysis of collec¬ 
tive oscillations 



Fig. 2 Illustration of the introduction of gaps in the zigzag 
arrangement in computer simulations. The gaps lead to longitudinal 
pairing cascades along the train of droplets.^ The average droplet 
spacing a in the configurations shown from top to bottom is 2.2R, 
2.AR, and 2.6R. 


Trajectories of 4 out of 40 droplets 



Fig. 3 Trajectories of 4 droplets in a train of 40 droplets obtained 
from simulations where the initial arrangement contains gaps. For 
the smaller spacing a = 2.2/? the longitudinal and transverse 
oscillations remain stable for some time until the longitudinal 
oscillations start growing. The transverse oscillations do not cross 
the channel centreline, and the trajectory in configuration space has 
a circular form. For the larger spacing a = 2.4/? the oscillations start 
to grow sooner and the droplets move transversely across the 
channel, thus breaking the initial pattern. 


of droplet i and no force is exerted due to the symmetry of the 
arrangement. The zigzag order is thus stable, and for a collec¬ 
tive oscillation to emerge the symmetry of the droplet arrange¬ 
ment has to be broken. In the following, we describe ways to 
excite oscillations in an array of droplets and analyse quan¬ 
titatively the collective modes observed in experiments and 
reproduced by computer simulations using MFC as described 
in section [3^ The results demonstrate the rich dynamics that 
emerge from hydrodynamic interactions in confined microflu¬ 
idic droplet ensembles. 


Arrays of droplets flowing in a straight microchannel self- 
organise into two parallel trains of droplets with alternating 
lateral positions.!^ For all neighbours j of a droplet i in this 
zigzag order, there exists a droplet / such that the positions 
relative to i satisfy • x = — • x and • y = • y. There¬ 

fore, the flow fields of droplets j and / cancel at the position 


4.1 Longitudinal oscillations 

In simulations, one way to perturb the symmetry of a droplet 
train is to vary the spacing by introducing ‘gaps’ between the 
incoming droplets in regular intervals. The resulting arrange¬ 
ment for the case of an alternating sequence of two differ¬ 
ent distances between subsequent droplets is illustrated in flg- 
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Fig. 4 Power spectra obtained from simulations of the oscillations 
triggered by the introduction of gaps. 


ure The droplet crystal now consists of pairs of laterally 
displaced droplets. For an isolated pair of spherical droplets, 
the hydrodynamic forces do not lead to relative particle mo¬ 
tion,^ owing to the time reversibility in the Stokes regime. 
In an ensemble, however, the pairs are also affected by the 
neighbouring pairs. In the perturbed zigzag arrangement, the 
pairs are separated by a larger distance than the droplets in the 
pair, and due to these gaps the configuration of the ensemble 
to the right of any droplet is different from the configuration 
on its left. Consequently, the flowing droplets can experience 
a net hydrodynamic force and undergo relative motion. We 
observe that the leading droplet of the pair moves faster than 
the trailing droplet. It separates and catches up with the trail¬ 
ing droplet of the pair ahead of the original one, thus forming 
a new pair. This process repeats with the new pairs and creates 
a longitudinal oscillation of the droplet distances in the train, 
where adjacent droplets are in anti-phase.^ Figureshows the 
motion pattern for this droplet arrangement in the co-moving 
frame of reference. The behaviour shows some similarity to 
the pairing cascades observed in finite droplet trains^SI, how¬ 
ever, in an infinite crystal the pairs cannot separate and keep 
interacting such that the oscillatory motion ensues. This is also 
reminiscent of the behaviour of colloidal particles driven by a 
constant external pulling force to move on a ring.®l The mo¬ 
tion patterns further reveal that in the initial stage, the trans¬ 
verse amplitude is small and the droplets stay on one or the 
other side of the channel. In this phase, the configuration 


Fig. 5 Space-time plots of simulation runs of a train of 40 droplets 
perturbed by gaps. The initial horizontal pattern indicates that the 
oscillations start out as a standing wave. After some time the initial 
pattern breaks, and waves begin to propagate along the droplet train. 


space trajectory of the droplets has a circular pattern. Over 
time, however, the transverse amplitude is growing and the 
droplets eventually cross the centre-line of the channel, as can 
be seen in figure[^for ajR = 2.4. 

The longitudinal and transverse power spectra of the droplet 
oscillations observed in simulations are shown in figure]^ For 
the smaller spacing a = 2.2/?, a clear signature of longitu¬ 
dinal oscillations is observed. The accompanying transverse 
oscillations are from droplets moving towards the centre of 
the channel and back without crossing the channel. For larger 
spacing a = 2.4/? the power spectra show more scattered fea¬ 
tures without a clear signature indicating the limited stability 
the gap modes. The longitudinal spectrum has a strong sig¬ 
nature CO oc k which is due to fiuctuations, but there are also 
distinct signals at k-a = 7l indicating a zigzag wave. 

The space-time plot of the droplet distance in figurej^shows 
for small times a pattern of horizontal stripes in the both the 
longitudinal 5x and transverse 5y displacements, which indi¬ 
cates that in the co-moving frame, the excited mode is initially 
a standing wave. For longer times, the pattern changes as a 
travelling wave seems to develop. These results suggest that 
the stability of the longitudinal pairing wave is limited, as we 
discuss in more detail below. 
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Fig. 6 a) Microscopy time series showing the droplet reorganisation 
at a 90° bend, b) Sketch indicating the relevant hydrodynamic 
interactions between the droplets in the co-moving frame. The 
transverse forces resulting from leading and trailing droplets in a 
zigzag configuration with one ‘mismatch’ are shown as grey and 
black arrows, respectively, c) Time series of a transverse droplet 
motion for a single ‘mismatch’ as described in b); three drops are 
marked and the net transverse forces, i.e., the sum of the grey and 
black arrows in b), are shown as black arrows. 


4.2 Transverse oscillations 

Another protocol to excite waves in the microfluidic crystal 
is to exchange the positions of a pair of droplets. Such ‘de¬ 
fects’ can be created experimentally by guiding the zigzag ar¬ 
rangement of droplets around a rectangular microfluidic bend, 
see flgure [^a). At the bend, the two droplets of a pair ex¬ 
change their longitudinal position as depicted. During this 
process, the passing droplet creates an accelerated flow at its 
trailing edge which prevents the following pair from exchang¬ 
ing positions. Under suitable conditions, every second pair 
undergoes a positional reordering such that the translational 
symmetry of the droplet train is broken. Hence, in the re¬ 
sulting droplet arrangement after the bend the droplets ex¬ 
perience a net hydrodynamic force. For certain droplet re¬ 
spectively channel dimensions, the bend allows to systemati¬ 
cally create such defects in the translational symmetry of the 
crystal, and if defects are created periodically a global oscil¬ 
lation patterns emerge. The collective oscillations are very 
stable and could be observed for channel lengths up to 10 cm, 
i.e. after travel distances which are four orders of magnitude 
larger than a typical droplet radius. The wavelength A in lon¬ 
gitudinal direction depends on both the droplet size and the 
droplet spacing. By a variation of these parameters, various 
initial wavelengths can be excited. Within a certain parame- 


0(^pOOC^OCPOCy^( ^CO^O 

V- v^'QjCO; 

Fig. 7 Experimentally observed travelling sine waves as generated 
by periodic ‘mismatches’ using the setup shown in figure [^having 
different droplet size {R6A jim, R61 jim, and R'^lOjim from 
top to bottom) and different wavelength (A = 8a, A = 6a, and 
A = 4a from top to bottom). 


ter range, we can achieve accurate control of the amplitude of 
the transverse oscillations by tuning the droplet radius R, since 
the amplitude of the transverse oscillations is equal to W — 2R, 
where W is the channel width. Figure |7] shows three examples 
of structures with variable wavelength and amplitude. They 
were obtained from zigzag structures with equal droplet spac¬ 
ing but with different droplet radii ^ 64 /im, R61 jim, 
and R ^ 70/im which represent the case of the largest, inter¬ 
mediate, and smallest wavelength producible. The first struc¬ 
ture presents 8 droplets, the second 6 droplets, and the last 
structure 4 droplets per wavelength. From the microscopy 
time series we extract the trajectories and from the trajecto¬ 
ries we could extract the transverse phonon spectra, see figure 

which reveal the presence of two peaks for the transverse 
oscillations. These peaks indicate that the excitations have 
distinctive wavelengths for both longitudinal and transverse 
modes. 

The same type of droplet wave also emerges in computer 
simulations when a triangle wave is used as initial condition, 
see figure Various wavelengths can be excited which are 
found to behave qualitatively the same. In the following, we 
analyse results obtained for wavelengths A = 5a and A = 6a, 
corresponding to five and six droplets within one wavelength, 
respectively. Figure[T^shows the trajectories of the droplets in 
the co-moving frame of reference. We find that both the lon¬ 
gitudinal and transverse coordinates oscillate around the equi¬ 
librium position. The configuration-space trajectory of each 
droplet describes a figure-eight pattern. While for A = 6a 
the figure-eight pattern is almost symmetric, it is clearly an¬ 
tisymmetric for A = 5 a. This is a consequence of an antisym¬ 
metric initial arrangement of droplets at the top and bottom 
walls where two neighbouring droplets are close to the top 
walls, while a single droplet is close to the bottom wall. The 
trajectories in figure \T0\ suggest that this asymmetry is main¬ 
tained in the oscillatory motion.^ The trajectories are reminis¬ 
cent of the bound-state motion pattern of an isolated pair of 
rigid discs in two-dimensional Poiseuille flow®. Quantita¬ 
tively, however, the state diagrams in Ref.® seem to predict 
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confinement i?/VF=0.30, spacing a/i?=3.29, X=8a 






Fig. 9 Snapshots of sine-like transverse waves observed in 
computer simulations starting from an initial triangle wave.^ The 
variable wavelength A in the configurations shown from top to 
bottom are 6 a, 5 a, and 4 a. 


confinement R/W=0.S2, spacing a/R=3.08, A=6a 



k-a k-a 



Fig. 8 Longitudinal and transverse power spectra of the droplet 
oscillations as calculated from experimental trajectories extracted 
from microscopy time series corresponding to the waves shown in 

Fig.|3 


cross-swapping trajectories for the relatively large displace¬ 
ments (|yi —y 2 \/y^ ~ 0.2 and 2a/W ^ 1.5) in the wave con¬ 
figurations considered here. We conclude that here the ensem¬ 
ble arrangement stabilises the oscillatory states, as already ob¬ 
served for the pairing waves as discussed in section [TT] for the 
cascade of droplet pairing. 

Figure m shows the longitudinal and transverse power 
spectra of the waves depicted in figure The transverse 
power spectrum from simulations shows a continuous sig¬ 
nature where the dependence of the frequency on the wave- 
vector has a sine-like shape. Such a dis persion relation is 
reminiscent of microfluidic phonons .ESEU xhe dispersion re¬ 
lation predicted by th e line arised far-field theory for confined 
microfluidic phonons, ®IE1 cf. section is plotted on top of 
the spectra and shows excellent quantitative agreement. It is 


Trajectories of 4 out of 90 droplets 

\ = f\n \ = ^n 



6x/a Sx/a 


Fig. 10 Trajectories of 4 droplets in a train of 90 droplets forming a 
sine-like wave as obtained from simulations for wavelengths X = 6a 
and X = 5a. The trajectories follow a figure eight-like pattern, 
which is symmetric for X = 6a and asymmetric for X = 5a. 


worthwhile to note that the calibration procedure described in 
section fixes all parameters such that no fitting is needed. 


The space-time plot of the droplet distance in figure [T^ con¬ 
firms that the wave is travelling in the flow direction along the 
droplet crystal, and is considerably more stable than the longi¬ 
tudinal waves discussed above. At long times, gaps start form¬ 
ing and the crystal breaks up into smaller sub-units. These 
results show that the excited transverse wave is an acoustic 
microfluidic phonon, and our experimental approach enables 
us to specifically excite such modes with a large amplitude. 
The power spectra also show a clear signature of longitudi¬ 
nal oscillations which is explained in detail in the following 
subsection. 


Here it is also possible to increase experimentally the am¬ 
plitude of the longitudinal waves by increasing the distance 
between the droplets. A droplet train with such ‘gaps’ is 
shown in figure along with the corresponding longitudi¬ 
nal and transverse power spectra. Compared to the sine-like 
waves without gaps, cf. figure the transverse modes here 
exhibit a broader signal around the main wavelength. More¬ 
over, the longitudinal signature extends over a whole range 
of wavelength, similar to the longitudinal modes observed in 
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confinement R/W=0.33. spacing a/R=2.22, X=6a 


Space-time plots of a train of 90 droplets, A=6a 



confinement R/W=0.SS, spacing a/R=2.22, A=5a 
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Space-time plots of a train of 90 droplets, A=5a 
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Fig. 11 Power spectra obtained from simulations of the sine-like 
waves with wavelength X = 6a (top) and X = 5a (bottom). The 
white lines are the dispersion relations {k) for transverse 
phonons, Eq. (fTSland the phenomenological dispersion relation 


ft)|| {k) from Eq. The theoretical dispersion relation (151 for 
longitudinal phonons is shown in grey. The dashed lines are the 
continuum approximations of Eq. (|24l. 


computer simulations, cf. figure]^ This indicates that a range 
of wavelengths have been excited longitudinally by increasing 
the droplet spacing, while the transverse sine-like wave is still 
predominant. From the optofluidic point of view, these hetero¬ 
geneous structures are the most interesting. When the droplet 
distance approaches the hydrodynamic screening length, more 
heterogeneous structures are observed in line with the results 


from computer simulations in section 4.1 


4.3 Coupling of longitudinal and transverse oscillations 

The experimental and numerical power spectra in figure and 
[^also show a clear signature of longitudinal oscillations, in 
particular for the simulation results. The dominant frequencies 
from experimental measurements agree quantitatively with the 
simulation results. However, the dispersion relation of the lon¬ 
gitudinal modes is qualitatively different from the prediction 
of the linearised theory. The frequency curve (k) for acous¬ 
tic phonons, c.f. Eq. ([Tg, is plotted as a grey line in figure 
[^and describes waves that propagate upstream, whereas the 
observed frequencies indicate a positive group velocity. The 
shape of the measured frequency curve resembles instead the 
shape of the transverse dispersion relation co±, yet the maxi- 


Fig. 12 Space-time plots of simulation runs of a train of 90 droplets 
forming a sine-like wave for wavelengths X = 6a and X = 5a. 


mum appears shifted towards the edge of the Brillouin zone 
and has a higher frequency than the maximum of co± (k ). 

The anomalous properties of the longitudinal modes are 
also apparent in the correlation strength 


C(a:||,/c_l) I- -- 

between longitudinal and transverse modes 


( 20 ) 


x\\{k^t) = ^ 5xj{t)exp 
f=i 

N 

yi{k,t) = Y, ^>'j(0exp 

;=i 




N 


( 21 ) 


The correlation strength for the transversely excited oscil¬ 
lations is shown in figure and indicates that longitudinal 
and transverse modes strongly correlate if and k^ have the 
same sign. This is in contrast to correlations that are expected 
for purely acoustic phonons, where the matching condition 
^11 = —ki_ is expected.®! Furthermore, by employing the fre¬ 
quencies for correlated wave-vectors in the observed disper¬ 


sion relations (figure 11), we find that the matching waves 
approximately obey mu =2co^. These observations indicate 
a strong coupling of the longitudinal modes to the transverse 


This journal is ©The Royal Society of Chemistry [year] 


Journal Name, 2010, [vol], 


























































^OMocP 


confinement R/W=0.3S, spacing a/R=3.01 



Fig. 13 (Top) Experimentally observed droplet pattern as generated 
by periodic ‘mismatches’ using the setup shown in figurej^with 
increased spacing or ‘gaps’ between some droplets. (Bottom) 
Corresponding longitudinal and transverse power spectra of the 
droplet oscillations as calculated from experimental trajectories 
extracted from microscopy time series. 


oscillation which is beyond a linearised far-field theory. This 
coupling is induced by the relatively large amplitude of the 
transverse waves which brings the droplets close to the wall 
where the imposed flow is not uniform but decays due to no¬ 
slip boundary conditions. ^3 The droplets are slowed down at 
both channel walls which leads to the flgure-eight trajectories 
seen in flgure[T^ In the co-moving frame, the droplets move 
backwards when they are close to the wall, and forward when 
they cross the channel. During one transverse cycle, two wall 
approaches take place such that 

( 0 || = 

in agreement with the matching condition found above. In¬ 
spection of the spatial wave patterns reveals that a full longitu¬ 
dinal wave extends over each crest or trough of the transverse 
wave, hence 

k|| = 2k^. 

Combining these conditions leads to the dispersion relation 





A=6a 


N=90 droplets 


A=5a 


Fig. 14 Correlation strength C{^,kj_) of phonon modes obtained 
from simulations of transversely excited oscillations with initial 
wavelengths X = 6a and X = 5a. The dashed line is the 
phenomenological relation ^ = 2k 




a R=2.20 


40 droplets 


a/R=2.A0 


Fig. 15 Correlation strength obtained from simulation 

data of longitudinal pairing waves triggered by gaps. 


We conclude that the lateral conflnement in connection with 
the considerable amplitude of the transverse excitations leads 
to a strong interaction of longitudinal and transverse waves in 
conflned microfluidic crystals. It is interesting to note that for 
the excitation mechanism studied here, this coupling does not 
seem to lead to an instability. A possible explanation is that 
in the dense droplet crystal, longitudinal oscillations cannot 
grow due to steric constraints between the droplets. 

4.4 Long-time behaviour and stability 


m||(k) =2(0 _l(V2). (22) 

This relation is also plotted in flgure and is in remarkable 
quantitative agreement with the signature of the longitudinal 
oscillations. 

The correlation strength for the longitudinal pairing waves 
discussed in section |4.1| is shown in flgure Here, the 
oscillations seem to be highly correlated in a region where 
^11 =k^ — 7l. However, this pattern was only found if the gaps 
between pairs were small, and since the pairing waves are con¬ 
siderably less stable, it is difficult to clearly identify the origin. 


While the transverse waves excited by positional exchange ap¬ 
pear to be relatively stable, the pairing waves excited by the 
introduction of gaps in the droplet train persist only for much 
shorter times. This is evident in the space-time plot of the 
longitudinal droplet distance in figure which extends the 
time scale of figure by 5 five times and shows that the ini¬ 
tial order disappears and other patterns emerge. One feature 
that is visible for a range of initial gap widths is the formation 
of subregions with small longitudinal droplet distances (dark 
red regions in the space-time plot). Inspection of configura¬ 
tion snapshots reveals that in these subregions, the droplets ar- 
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Space-time plots of a train of 40 droplets, a/R=2.AQ 
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Using this continuum approximation to derive the equations of 
motion and the dispersion relation as in section]^ we obtain 

(JL)||(^) = —4^(2^ coth(7;rj3)csch^(7;rj3) 

co±{k) = 2kaB [1 +cosh( 7 ;rj 8 )]^csch^( 7 ;rj 8 ). 

These relations are plotted as dashed lines on the power spec¬ 
tra in figure [TT] The longitudinal pairing waves clearly show 
a dispersion branch that is linear in k, and the spectra of the 
transverse waves also have a signal that agrees with a linear 
dispersion relation. We find that this part of the spectrum 
becomes more pronounced at longer simulation times, which 
supports our hypothesis that the droplet train becomes disor¬ 
dered due to fluctuations. Overall our analysis sheds some 
light on the long-time behaviour of oscillations in microfiu- 
idic crystals, however, more data will be needed to study the 
onset and growth of instabilities further. 

5 Summary and discussion 


Fig. 16 Space-time plots of simulation runs of a train of 40 droplets 
initially perturbed by gaps with average spacing between the 
droplets a/R = 2.2 and alR = 2 A, respectively. 


range in a dense zigzag order without oscillations. We hypoth¬ 
esise that the formation of these ordered regions as observed 
in simulations is similar to fiow-induced crystallisation.^SThe 
‘frozen’ parts of the droplet train propagate with a velocity 
that depends on the size of the subregion, and due to differ¬ 
ent velocities some zigzag clusters catch up and merge with 
others. Simulation trajectories reveal that the droplets crys¬ 
tallise at the leading edge and melt at the trailing edge of the 
subregions. In between the ordered regions, the droplet train 
is disordered. In some instances, we observe a phonon-like 
transverse wave that develops in these regions.^ One simula¬ 
tion snapshot of a droplet train that simultaneously exhibits an 
ordered zigzag region and a transverse phonon wave is shown 
in figure [T^ 

In the simulations, the long-time behaviour is significantly 
affected by thermal fluctuations, which tend to lead to a disor¬ 
dered droplet ensemble where the crystal structure disappears. 
When the droplets move away from their regular crystal posi¬ 
tions, they can be regarded as a continuous ensemble, and the 
finite-differences in the equati on of motion are to be replaced 
by a continuum approximation®!!! for small a 


R e ^ . d5y 

oyn^j - oyn-j - 2ja-^. 


(23) 


Excitation mechanisms for collective waves in microfiuidic 
crystals have been investigated. We have demonstrated that 
both longitudinal and transverse waves can be systematically 
excited by creating specific defect patterns. Experimental re¬ 
sults were confirmed by computer simulations, and our results 
reveal instabilities and mode coupling that originate from the 
underlying hydrodynamics. 

The excited longitudinal modes show cascades of pairs 
of laterally displaced droplets. Due to the pairing instabil¬ 
ity, the pairing cascade is rather short-lived, and over time 
other modes are observed. Since the constraints prevent the 
dense crystal from becoming completely disordered, some 
parts ‘freeze’ into a zigzag arrangement while others exhibit 
transverse oscillation. This is a possible indication of fiow- 
induced crystallisation, and it will be interesting to further 
investigate the dynamic formation of the inhomogeneous re¬ 
gions. 

The transverse waves show the dispersion relation of a mi¬ 
crofiuidic phonon. Comparison with the analytical prediction 
demonstrates that a linearised far-field theory works well even 
in a dense droplet crystal. A possible reason is the relatively 
strong lateral confinement, which may screen the higher-order 
reflections of the dipolar interactions. The power spectra of 
the oscillations exhibit a correlation that arises from a cou¬ 
pling of longitudinal to transverse modes. This coupling is 
induced by the boundary conditions at the confining channel 
walls. At large amplitudes, the inhomogeneity of the imposed 
flow affects the dynamics of the droplets and leads to novel 
collective modes. Our results thus shed light on the mech¬ 
anisms underlying non-linear mode coupling in microfiuidic 
crystals. 
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Fig. 17 Simulation snapshot illustrating the formation of ordered and unordered subregions within a train of 40 droplets. The right part of the 
train has ‘frozen’ into a stable zigzag configuration, while the left part exhibits a travelling sine-wave.^ 


Experimentally, the excited waves are highly stable and do 
not undergo instabilities. This raises the question wh ether 
some of the instabilities observed in microfluidic crystals^ini 
are suppressed in the dense droplet ensembles studied here. 
Large displacements are inhibited by steric interactions be¬ 
tween the droplets which may have a stabilising effect, e.g., 
on the pairing cascades. A detailed study of the impact of ge¬ 
ometric constraints on the stability of the collective modes is a 
topic for future research. In the simulations, the smaller Peclet 
number leads to a more pronounced influence of thermal fluc¬ 
tuations. This opens up the possibility for fluctuation-induced 
instabilities and our simulation approach may thus be used to 
investigate instabilities in microfluidic devices. 

The good agreement of the experimental and numerical re¬ 
sults with the linearised far-fleld prediction suggests that the 
dynamics of dense and conflned microfluidic droplets is acces¬ 
sible theoretically and leads to novel insights into the origin of 
instabilities and mode coupling effects. Furthermore, it seems 
promising to apply the experimental and numerical techniques 
to other microfluidic systems, such as dense droplet systems 
in 2D, crowded particle systems, or self-propelled particles. 
Finally, the experimental techniques may be used to control 
particle flows in microfluidic applications such as flow cytom¬ 
etry or high-throughput assays using microchips. 
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